function [rkstat, pval_Wald, pval_F]=rank_eig_multiasset(R, F,assetclasses)
%{
    This function implements the rank test (eigenvalue approach)
    
    INPUT 
    R       : n by T, a matrix of asset returns
    F       : k by T, a matrix of risk factors
    assetckass: n by 1, integers 1 through nassets indicating asset classes
    
    OUTPUT  
    rkstat  : rank test statistic


    NOTATION: |\theta*B-A|=0, see Kleibergen (2009)
%}

    nassets = max(assetclasses);

    [n, T]=size(R);                       % get n, T
    k=size(F,1);                          % get k
    F_bar=F-mean(F,2)*ones(1,T);          % derive F_bar
    B=inv(F_bar*F_bar'/T);
    
    Rs = zeros(1,T);
    
    for j=1:nassets
    
        assetsj = assetclasses == j;
        nj = sum(assetsj);
        if nj > 1
            Rj = R(assetsj,:);

            ln=ones(nj-1,1);                       % a vector of ones
            Rsj=Rj(1:(nj-1),:)-kron(ln,Rj(nj,:));      % delete the n_th asset    
            Rs = vertcat(Rs,Rsj);
        end
        
    end    
    
    Rs = Rs(2:end,:);
    
    %writematrix(Rs','./MatlabInput/test.csv');
    
    F_exp=[F;ones(1,T)];
    B_temp=Rs*F_exp'/(F_exp*F_exp');
    Sigma_tt=1/T*(Rs-B_temp*F_exp)*(Rs-B_temp*F_exp)';
    Bhat=B_temp(:,1:end-1);
    
    A=Bhat'*inv(Sigma_tt)*Bhat;           % matrix A
    
    % rank test
    [~,D]=eig(A,B);                       % V:eigenvectors; D:eigenvalues  
    rkstat=T*min(diag(D));
    %n-nassets+1-k
    pval_Wald=1-chi2cdf(rkstat,n-nassets+1-k);
    cRank=k-1;
    dg1=(n-nassets-cRank)*(k-cRank);	
    dg2=(T-k-1)-dg1+1;
    stat_F=rkstat/T*dg2/dg1;
    pval_F=1-fcdf(stat_F,dg1,dg2);
end 